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The energy of atmospheric neutrinos detected by MACRO was estimated using multiple coulomb scattering 
of upward throughgoing muons. This analysis allows a test of atmospheric neutrino oscillations, relying on the 
distortion of the muon energy distribution. These results have been combined with those coming from the upward 
throughgoing muon angular distribution only. Both analyses are independent of the neutrino flux normalization 
and provide strong evidence, above the 4a" level, in favour of neutrino oscillations. 
PACS: 14.60Lm,14.60.Pq, 25.30.Mr 



1. Introduction 



The results obtained by experiments looking 
for neutrinos coming from natural or artificial far 
sources gave the first indication of physics beyond 
the standard model. Neutrino flavor changing is 
in fact the most straightforward explanation for 
electron and muon neutrino disappearance and 
for the evidence of non electron neutrinos from 
the sun. In this scenario, the study of atmo- 
spheric neutrinos plays an important role to im- 
prove our understanding of the neutrino oscilla- 
tion mechanism. 

The MACRO detector studied [TEin) three cat- 
egories of events as shown in Fig. 1: (1) upward 
throughgoing muons, (2) upward semicontained 
muons, and (3-1-4) upward stopping muons plus 
downward semicontained muons. The (3-f 4) cat- 
egories cannot be separated experimentally, be- 
cause the events have only one time measurement. 
However, upward stopping muons and downward 
semicontained muons have similar parent neu- 
trino energy. 

An atmospheric muon neutrino deficit was 
found in all of the three categories and an angu- 
lar distribution distortion was found in category 
(1). Such results are explained by the neutrino 
oscillation hypothesis with parameters Am^ = 
2.5 X 10~^eV^ and sin^20=l, in good agree- 
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Figure 1. Crp&§ section of the MACRO detector 
and topologies of events induced by neutrino in- 
teractions. 

ment with the Super-Kamiokande results 
A detailed study of the upward throughgoing 
muon angular distribution by MACRO ^ and 
by Super-Kamiokande [S] allowed the exclusion 
at the 99%C.L. of the muon neutrino into sterile 
neutrino oscillation mechanism, compared to the 
muon neutrino into the tau neutrino. 

Considering a two flavor neutrino oscillation, 
the probability for a Vfj_ to oscillate to is given 
by: 

P,^^,^ = sin2 26 sin' j (1) 

where L^{km) is the distance between the neu- 
trino production and interaction points and 
i?i/(GeV) is the neutrino energy. The different 
categories of events quoted above have median 
energy from 4 GeV to 50 GeV, which provides ev- 
idence of the dependence on neutrino energies, as 
required by the oscillation hypothesis. In this pa- 
per we address for the first time the estimate of 
the upward throughgoing muon energy by using 
multiple scattering. The implementation of this 
method to the MACRO data is discussed in Sec- 
tion 2. Two different analysis were performed. 
The first one, using the tracking system in digital 
mode, is described in Section 3. Section 4 shows 
the results obtained with the streamer tube in 



drift mode, where electronic readout of the hit 
time was used to improve the position resolution 
and thereby the range of muon energies that could 
be estimated. 

2. The muon energy estimate 

The MACRO detector was extensively de- 
scribed in 0. The detector consisted of a lower 
half (the lower detector) and an upper half (the 
"Attico"). The lower detector had 10 horizon- 
tal planes of streamer tubes separated by either 
crushed rock absorbers (60g/cm^) or (at its top 
and bottom) a plane of scintillator counters. The 
Attico had 4 horizontal planes of streamer tubes 
separated into two groups with a plane of scintil- 
lator counters between them. 

Upward throughgoing muons are mainly pro- 
duced in neutrino deep inelastic scattering (DIS) 
in the rock below the detector. The recoil hadrons 
are lost and the muon energy is degraded in the 
propagation to the detector. Nevertheless, Monte 
Carlo simulations show a linear relation between 
the parent neutrino energy Ei, and the muon 
residual energy at detector level. 

The momentum resolution obtainable from 
multiple coulomb scattering (MCS) measure- 
ments is the result of two different contributions: 
the number of sampling planes A'' and the space 
resolution of the detector tJo- In MACRO the 
number of tracking planes interleaved with rock 
absorbers gives A^ = 8. The resolution of the 
streamer tube system used in digital mode is 
ctq — 1 cm and was improved to ctq ~0.3 cm by us- 
ing the drift time. The other six horizontal track- 
ing planes are separated by a negligible amount 
of material, and do not contribute to the MCS 
measurements. 

In a tracking detector with equispaced track- 
ing planes, separated by a given absorber and 
neglecting the energy losses, the characteristic 
momentum scale can be written as 
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1><PmcSj the main limitation to the momentum 
reconstruction comes from the number of sam- 
pling planes while for p>Paics the space reso- 
lution dominates. Under the conditions quoted 
above the relative momentum error |8j can be 
written as: 
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P a/2N Vpmcs, 

for p>>Pmc'S- As shown in [5] the 1/p distri- 
bution approaches a Gaussian only for N>30; 
therefore the momentum error for MACRO is 
not expected to have a Gaussian behaviour. The 
eight horizontal streamer tube planes interleaved 
with rock absorbers (planes 2 to 9, from the bot- 
tom) are equispaced and have the same space 
resolution, hence p])^^g^°=2.2 GeV/c. Since 
about 90% of upward throughgoing muons have 
p > 2.2 GeV/c, the intrinsic chamber resolution 
dominates. The MCS-based momentum esti- 
mates scale therefore as the square of the space 
resolution CTq- 

The r.m.s. of the lateral displacement of a rel- 
ativistic muon crossing a layer of material with 
depth X is proportional to the inverse of muon 
momentum |^ : 
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Xo, 



where A is the distance between tracking planes 
and Xq is the material radiation length. For 



In MACRO, X~25Ao/cose, giving on the verti- 
cal a^^^f ~10cm/£;^(GeV). Therefore, a satu- 
ration point above 10 GeV requires a space reso- 
lution better than 1 cm. An improvement in the 
space resolution increases the maximum energy 
value (saturation point) above which the MCS 
method is not effective. 

Consequently we decided to improve the 
streamer tube resolution: in the second analysis, 
described in Section 4, we took advantage of the 
QTP-TDC electronics, developed for magnetic 
monopole searches [7]. In this case, we obtained 
an improved space resolution, i7^~0.3 cm|ll|. 
which allowed us to estimate the muon energy 
up to -40 GeV. 



For both analyses, we used the whole sample 
of upward throughgoing muon events collected 
with the complete apparatus in a period of data 
taking equivalent to 5.5 years of live time. We 
studied upward muon events selected by the time- 
of-flight measured by planes of scintillators com- 
bined with the standard MACRO tracking algo- 
rithm. The original upward throughgoing muon 
data set has been described previously To 
make a comparison between real data and ex- 
pectation we performed a Monte Carlo simula- 
tion using the Bartol neutrino flux ^2] and the 
GRV94 DIS parton distributions ^] for deep in- 
elastic scattering. For low energy neutrino inter- 
actions we used the cross sections given in |14) . 
The energy loss for muons propagating through 
rock is taken from Lohmann et al. jl5j while 
the muon simulation inside the detector was per- 
formed with GMACRO (the GEANT 3.21 based 
detector simulation). For the second analysis, 
where the streamer system is used in drift mode, 
another simulation chain has been implemented. 
In this case neutrino interactions are randomly 
distributed in a rock semi-sphere below the de- 
tector. Muons are then transported to the de- 
tector using the FLUKA99 package 6^. A Monte 
Carlo statistics corresponding to a live time of 
~2700 years was produced (500 equivalent exper- 
iments). This simulation was compared with that 
used in 1.3 , obtaining a satisfactory result. 

3. The first analysis: streamer tubes in 
digital mode 

Because of MCS in the lower detector, with 
this analysis we expect a measurable deflection 
between the incoming and the outgoing track di- 
rections for muons with energies smaller than 
~10GeV. This corresponds to an effective arm- 
lever of ~4 m between the lower detector and the 
Attica^. 

The eight lowest streamer tube planes were 
used, through a track refit, to estimate the di- 
rection of the incoming muon. The five upper 
streamer tube planes in the lower detector and 
the four Attico planes were used to estimate the 
direction of the outgoing muon. The distance r^, 
between the intercepts of the two tracks in the 



z — plane, and the difference A$ between the 
two slopes depend on the muon energy E^. We 
divided the upward throughgoing muons in three 
sub samples, according to the values of and 
A$: sample L= Low energy, if r^ >3cm and 
A$ > 0.3°; H= High energy if r^ <3cni and 
A$ < 0.3°. The remaining events were classi- 
fied as M— Medium Energy. These cuts were opti- 
mized using two large samples of real atmospheric 
muons: i) downward throughgoing (average en- 
ergy of ^ 300 GeV .16 ), and ii) downward going 
stopping muons (average energy of ~1 GeV). 

The aforementioned cuts were applied on the 
upward throughgoing muons (crossing the whole 
apparatus, lower detector + Attico), both real and 
simulated. They were selected from the same 
analysis chain ^ and had the same data for- 
mat. From the simulation we expect 430 events: 
178 of Low energy, with <£;^>=llGeV, 59 of 
Medium energy, with <i?;^>=33GeV and 193 of 
High energy, with <E,y>—72GeV. From the real 
data, 316 events were selected: 101 of Low, 51 of 
Medium and 164 of High energy. 

The L and H events were further divided ac- 
cording to their zenith angle 0: events from the 
vertical direction —1. < cos0 < —0.8, and events 
with COS0 > —0.8. For each event topology the 
number of detected and expected events were de- 
termined and ordered with decreasing value of the 
<Ly>/ <E,y> (<Ly>= average value of neutrino 
path length). The relative systematic uncertainty 
on each of the five <L,y>/<E^> values is 12.5%, 
coming from the uncertainties on the neutrino 
spectrum and angular shape (discussed in Sect. 
4.3), from the detector related effects and analy- 
sis cuts. This includes the number of planes used 
for the refit, the definition of the cuts, the fluc- 
tuations in the streamer tube and scintillator ef- 
ficiencies, and detector acceptance uncertainties. 

We evaluated for the hypothesis of no neu- 
trino oscillations using the five <L^>/<E,y> 
bins, plus the additional point from the analysis 
of semi contained upward going events (lU) [5] 
(upgoing neutrinos with <i?i,>~4 GeV). We 
found that the distribution agrees with the no os- 
cillation hypothesis with a probability lower than 
2%). For two-fiavor oscillations with parameters 
Am^ = 2.5 X 10~^ eV^ and sin^ 26* = 1 we get a 



probability of 45%. 



4. The second analysis: streamer tubes in 
drift mode 

The performance of the streamer tube system, 
read out in drift mode, is described in Here 
an absolute muon energy calibration was per- 
formed at the PS-T9 and SPS-X7 beams, where 
a slice of the MACRO detector was reproduced. 
The MCS information was handled by a neural 
network, based on JETNET 3.0[ni, calibrated 
with the muon beams quoted above. The neu- 
ral network(NN) output obtained using test beam 
data, was compared with that expected from the 
Monte Carlo simulation, obtaining a satisfactory 
agreement^!]. The application of such analy- 
sis to the MACRO data resulted in an improve- 
ment of the space resolution from Cx— 1 cm to 
ax—O-i cm. 

Figure 121 shows the Monte Carlo prediction for 
NN downward throughgoing muons and down- 
ward going stopping muons compared with ex- 
perimental data. A nice agreement between sim- 
ulation and real data is found for both categories, 
representing a large energy interval. 

The result of the neural network output energy 
calibration shows that NN output is almost linear 
with logiQ{Efj_), increasing with the muon energy 
up to Efj,~40 GeV, where a saturation effect oc- 
curs. With respect to the approach used in[TT] 
few details of the neural network were optimized: 
the neural network training and the energy cali- 
bration was performed separately for events with 
hits in the upper part of the detector (j4ttico). 

Although smeared by the energy carried away 
by hadrons and by energy loss in the rock, the de- 
tected neutrino induced muons still carry on av- 
erage ~40% of the original neutrino energy. By 
using the full Monte Carlo simulation quoted in 
Sect. 2, we calibrated the NN output as a func- 
tion of log]^Q(£'ix): the calibrated NN output is lin- 
ear up to log2o(£'^(GeV))~2.15. We may write: 

SE 

6 logio(-E^) = logio(e) ^In^;,, = logio(e)— ^, (6) 

where d\ogiQ{Ei,) is the difference between the 
logio of the real and of the reconstructed neutrino 




Figure 2°. Neural network(U£L)a(}^mtiJ^irA9wn- 

throughgoing muons: Real data (black squares) 
and Monte Carlo expectations (continuous line). 
NN output for downward going stopping muons: 
real data(empty circle) and Monte Carlo expecta- 
tions(dotted line). 

energy. Taking into account that the NN out- 
put was calibrated as a function of \ogiQ{E^), the 
energy resolution 5E^/Ey was obtained by plot- 
ting the quantity 5 logiQ(£'j.)/logiQ(e). In Fig. |31 
we show the resolution that could be obtained 
with an ideal muon energy resolution(dotted 
line) and that obtained with the present analy- 
sis (continuous line). The precision of the neu- 
trino energy estimate obtained with an ideal 
muon energy resolution detector is SEi,/E^~70%, 
while with the present method a resolution of 
SE^/Et,~150% is obtained. The asymmetry 
present in both curves comes from neutrino inter- 
actions occurring far from the detector, for which 
a large fraction of the muon energy is lost during 
the transport. 

4.1. Data Selection 

We used the 783 upward throughgoing muon 
data set collected in the full detector run started 
in 1994. Table 1 describes further event selection 
to arrive at the sample for MCS analysis. We re- 
quired a single track in both the wire and the strip 
view. We selected hits belonging to the track and 
made of a single fired tube, to associate unam- 
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Figure 3. NeAitrino en<j£g?^^esaja^i^™t^/irig,q(^ he 

obtained with an ideal residual niMon energy reso- 
lution (dotted line) and with the MACRO energy 
estimate (continuous line). 



biguously its QTP-TDC time information. This 
cut is effective for muon tracks with large zenith 
angles (0 > 30°), while it is quite loose around 
the vertical: we restricted the present analysis to 
events with G<60°. Comparisons were performed 
between simulated and experimental downward 
going muon data, to ensure that the selection ef- 
ficiency, in the used angular window, is the same 
in the Monte Carlo and in the real data. Spu- 
rious background hits have been avoided by re- 
quiring a time window of 2 jis around the trigger 
time. Finally, we selected events with at least 
four streamer tube planes with valid QTP-TDC 
data. We fitted the drift circles using the same 
tracking developed to analyse test beam muons. 
A minimum path length of 200 cm in the lower 
detector is required for tracks hitting the Attico 
and 400 cm for tracks not hitting it. These ge- 
ometric requirement ensure a minimum depth of 
material where muons may experience a measur- 
able amount of multiple scattering and a lever 
arm long enough for a comfortable tracking. Af- 
ter the selection cuts 300 events survived, giving 
an overall efficiency of 38.3%. 



Cuts 


Number of 
events 


Total number 


783 


of upward 




throughgoing muons 




Single track in the 


695 


wire and strip views 




>4 planes with 


347 


valid TDC hits 




Track length cut 


314 


e<60° 


300 



Table 1 

Data selection cuts and selection efficiencies. 



4.2. Qualitative tests 

We used the information provided by the neu- 
ral network to separate the neutrino events into 
four energy subsamples, as shown in Table 2. The 
same selection was applied to simulated events. 

Figure 4 shows the zenith angle distributions 



Sample 


Energy cuts 
(GeV) 


Median energy 
(GeV) 


Low 


£^;;'='=<30 


13 


Medium-Low 


30<.B^«=<80 


36 


Medium-High 


80<£;^'='=<130 


88 


High 


i;^^'=>130 


146 



Table 2 

Subsamples selected according with the reconstructed 
neutrino energy. 



of the upward throughgoing muons in the four 
energy subsamples, compared with the expecta- 
tions of the Monte Carlo simulation, assuming 
no-oscillations(dotted line) and oscillations with 
parameters Am^= 2.5x10"^ eV'^ and sin^2^ ~1 
(wall boxes). The Monte Carlo including the os- 
cillation hypothesis, with the parameters quoted 
above, reasonably reproduces the real data in 
each subsample. We point out the strong differ- 



<IJ 35 -- 



1 10 ~- 



<U 40 ~= 




the zenith angle (^^%&four energy ranges, 
points are the real data, dotted line is the Monte 
Carlo simulation, assuming no oscillation, and 
dotted boxes are the Monte Carlo expectation with 
Am^ =2.5x 10"^ eV^ ands\v?26 = I, including 
a 11% error. 

ence between the oscillations/no-oscillations hy- 
potheses at low energies, while such difference is 
reduced by increasing the reconstructed neutrino 
energy. 

Finally, we used information on the ratio 
L^/E^. The distance L^, travelled by the neu- 
trinos from production to the interaction points, 
was measured by MACRO relying on the muon 
zenith angle determination, with a precision 
AZ/y/Z/^ ^ 3%. The resolution on the ratio 
LyjEy is therefore fully dominated by the uncer- 
tainty in the neutrino energy estimate, giving a 
relative error of ~ 150%. The ratio DATA/Monte 
Carlo as a function of \o%-^Q(Ly j Ey), is plot- 
ted in Fig. |S1 The black circles are the real 
data over the Monte Carlo predictions, assuming 
no oscillations, the shaded regions represent the 
Monte Carlo predictions with oscillations, Am^ 
= 2.5xlO-3eV2 and sin^ 26^=1, divided by the 
Monte Carlo predictions with no oscillation. 

The log]^o(^) distribution of the neutrinos de- 
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Figure 5. Ratio (Data/MC) as a Junction of 
the estimated Ly/Ey for the MACRO upward 
throughgoing muon sample. The black circles 
are the real data over the MC(no oscillation), 
the solid line is the MC assuming Am^ = 2.5x 
10^^ eV"^ and sin^ 29=1 over the MC with no os- 
cillation. The shaded region represents the MC 
uncertainties. The last point (empty circle) is ob- 
tained from semicontained upward going muons. 

tected by MACRO spans from 0.8 to 3.5. The 
left-pointing arrow at low log]^Q(^) represents 
the effect of the neural network saturation. Since 
the maximum energy reconstructed by the NN is 
ii^i/~140GeV and since the minimum Ly for the 
present analysis is LJJ'™— 6500 km, the minimum 
estimated value is log]^Q(^)"""~1.7. As far as 
the right-pointing arrow at high ^ is concerned, 
the reconstruction of neutrino energy based on 
the residual muon energy, does not allow one to 
reconstruct the ratio L/E beyond a given limit. 
This is due to far neutrino interactions, originated 
by high energy muons, which lost a large fraction 
of its energy. An ideal detector with ideal muon 
residual energy resolution, requiring E^ <2 GeV, 
would select an average logiQ(^)=3.2. Finally, 
the neural network resolution in muon energy es- 
timate results in a maximum value log]^Q(-g^)=3. 
The dashed line at 1 is the expectation without 
oscillations. The last point (empty circle) is ob- 
tained from semicontained upward going muon 
rate|2]. It is not used for the evaluation of the 
oscillation probabilities (Sect. 4.4) nor in the al- 
lowed region plot (Fig. 7). Good agreement is 
found with the oscillations expected with the pa- 



rameters quoted above. 



4.3. Experimental results 

To quantify the contribution of the multiple 
scattering measurement in stand-alone mode, we 
performed a blind analysis using the Monte Carlo 
data, looking for the variable, based on the en- 
ergy estimate, showing the maximum sensitivity 
to separate the oscillation from the no-oscillation 
hypothesis. We found that the best performance 
is given by the ratio: 

R ^ Nion,lNMgh (7) 

where Nio^ and Nhigh are the number of events 
with £;7'=<30GeV and .B^"==>130 GeV respec- 
tively. 

We considered systematic uncertainties due to 
the neutrino flux calculation and neutrino cross 
section. Due to the large uncertainty in the abso- 
lute flux we use in this paper only 

the angular distribution and the ratio between 
different event categories selected according to 
the reconstructed energy. Nevertheless, the the- 
oretical uncertainty of the existing neutrino flux, 
based on the old CR spectrum 12 , has to be ac- 
counted for. We varied the input primary cos- 
mic ray spectral index 7 in our simulation by 
A7=±0.05, obtaining a theoretical error on R, 
(^)/z..=±13%. 

Another source of systematics comes from the 
neutrino cross section. We looked at R vary- 
ing the Monte Carlo input cross sections. The 
most important contribution to the error comes 
from the "Low" category, since the cross section 
at low neutrino energy is more uncertain. By 
comparing the cross section computed under sev- 
eral hypotheses, varying for instance the struc- 
ture function (521,1231) in the deep inelastic scat- 
tering, including or neglecting the contribution 
of resonant scattering, we found {^^) ^—^ %■ 
We estimated a total theoretical error (^w^\, = 

\ H / theo 

\/(Tr)k. + (Tr)'=16^- The systematic er- 
ror on R reconstruction, evaluated in 6%, includes 
the uncertainties on absorber density, drift veloc- 
ity, streamer tube efficiency and detector accep- 
tance. 
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Figure 6. Low energy events over high energy 
events as a function of Am^ .• the area between 
the solid lines includes a 17% systematics (the er- 
ror is not gaussian, see text). The hatched band 
represents the real data. 

Figure El shows the ratio R as a function 
of Am^, assuming maximal mixing, for the 
Monte Carlo simulation: the area between the 
two solid curves represents the 17% system- 
atic error. The Monte Carlo prediction in case 
of no oscillations is i?'=1.5±0.25(theo.-|-sys.), 
while for Am2=2.5x lO^^eT^^ and sin^ 26* =1, 
R=1.00±0.17(theo.-t-sys.). As pointed out inO, 
the ratio does not have a Gaussian distribution. 
The errors on the ratio are therefore also not 
Gaussian: they are reported just to give a crude 
estimates of the significance. The experimental 
ratio is i?'=^P=0.85±0.16(stat.). The one sided 
probability of measuring a value smaller than the 
measured one, was computed by using two differ- 
ent methods. In the first one we let the Monte 
Carlo simulation (no oscillations) to fluctuate ac- 
cording to statistical and systematic errors of the 
considered ratio. We then evaluated the fraction 
of events giving a value smaller than the mea- 
sured one. In a more pessimistic view, rejecting 
the hypothesis that the lower number of events 
detected by MACRO with respect to the Monte 
Carlo expectation is due to oscillation effects, the 
total number of Monte Carlo events (sum of the 



numerator and the denominator in the ratios), 
were normaUzcd to the experimentally measured 
sum and then they were allowed to fluctuate. 

The corresponding one sided probability of 
measuring a value smaller than R'^^p according 
to method 1 (2) is 0.75%(1.9%) corresponding to 
2.4cr(2.1 fj). Finally, we combined the informa- 
tion coming from the energy estimate and from 
the angular distribution. Wc considered the ratio: 



A = N^ert/Nfior, 



(8) 



where Ny^rt is the number of upward through- 
going muon events with cos0<-O.7 and Nhor is 
the number of events with cos0>-O.4, as dis- 
cussed in [3]. The probability of measuring 
a value A' lower than A^^^ according to the 
method 1 (2) is P(A'<A'==^p)=0.001%(0.01%), 
corresponding to 4.3 a{3.7a). It is worth not- 
ing that, combining the two independent prob- 
abilities on the ratios R and A, the probability 
that a fluctuation of the expected values (no os- 
cillations) generates the experimental results is 
P{^'-i?'<^'=^P-i?^^P)=1.3x 10-^(3. 2x 10-^), cor- 
responding to 4.7 o-(4 a). 

The 90% confidence level allowed regions in the 
oscillation parameter space have been computed 
according to the prescription given in j24j . Fig- 
ureHshows the 90% C.L. for the ratio R, the an- 
gular distribution ^ and for their combination. 

5. Conclusions 

Muon multiple coulomb scattering has been 
used to estimate the energy of the neutrino in- 
duced upward throughgoing muons detected by 
MACRO. The different tests performed on the 
data, using the estimated energy (separation in 
subsamples with different energies, L,y/E,y esti- 
mates, etc.) give a consistent picture, all of 
them supporting the neutrino oscillation hypoth- 
esis with parameters Am^ = 2.5xl0^"^ey^ and 
sin^ 20~1. To quantify such effect, the ratio 
R = Niow I Nhigh was used, in stand-alone mode 
and in combination with the angular distribution, 
A — Nyert/Nhor- Both of them are independent 
of the neutrino absolute flux. The significance of 
the MACRO observation of the neutrino oscilla- 
tions is above 4cr. 



1: Angular distribution 
2: Energy(Low/l-ligli) 
1+2 




Figure 7. 90% C.L. allowed region in the (Ss^<) 
siv?29) plane for v^^Vr oscillations, obtained 
with different data samples. 
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